function [msgOut,MtsubOut,tsOut] = writeNeuMom(filena,OFP,ts,Mos,strData,Mtsub)
%writeNeuMom - Write one or more .neu/.mom files for CATS/Hector 
%
%   [msgOut,MtsubOut,tsOut] = writeNeuMom(FILENA,OFP,tsIn,Mos,StrData,Mtsub)
%
% This fuction generates ASCII file(s) with GNSS data suitable for further 
% processing with CATS or Hector.
% Options on input arguments:
%
% FILENA is the name of the generated ASCII file. Options:
% - If FILENA is a char variable, it is the name of the output filename.
%   The saved file has extension '.neu' or '.mom' according to the choice 
%   about further data processing ('.neu' for CATS, '.mom' for Hector). 
%   The output file carries the three components in the CATS case, whereas
%   the three files '(FILENA)_0.mom', '(FILENA)_1.mom', '(FILENA)_2.mom'
%   are generated in the Hector case for the East, North and Up direction
%   respectively. 
%   If the time subdivision is carried out (see below), the filenames are 
%   '(FILENA)_1.neu', ... '(Filena)_M.neu' (CATS case) or '(FILENA)_1_0.mom',
%   '(FILENA)_1_1.mom', '(FILENA)_1_2.mom', ..., '(FILENA)_M_0.mom',
%   '(FILENA)_M_1.mom', '(FILENA)_M_2.mom', where M is the subdivision 
%   number;
% - If FILENA is undefined or empty, it is FILENA = statName, where statName
%   is the value of the field statName of the object tsDataIn.
%
% OFP is the option for further processing. Options:
% - if OFP is undefined, empty or false, the further processing is assumed 
%   to be carried out with CATS;
% - in all the other cases the further processing is assumed to be carried 
%   out with Hector. In particular, 
%   -- if OFP is 1 or the vector [0 1 2], three ASCII files are generated 
%      for each time segment (convention: 0 East, 1 North, 2 Vertical);
%   -- if OFP is the vector [0 1], two files are generated for each time
%      segment (horizontal components);
%   -- if OFP is the scalar 2, a file (vertical component) is generated
%      for each time segment;
%   -- in each other case, under the condition that OFP is true, it is as
%      if OFP is 1.
%
% tsIn is the input tsData object. Options:
% - If tsIn is a tsData object, it is used to take all the data about
%   the managed GNSS station.
% - If tsIn is undefined, empty or is not a valid tsData object, the
%   name of a MATLAB .mat file is interactively managed. Such a file should  
%   carry one or more variables. If a single variable is carried, it must 
%   be a tsData object. If two or more variables are carried, one and only 
%   one of them must be a 'tsData' object.
%
% Mos is the (optional) offset matrix. Options:
% - If Mos is a No-by-2 matrix, Mos(k,1) is the k-th offset date, expressed
%   in decimal form, and Mos(k,2) indicates the corresponding component(s) 
%   according to the "binary notation" used by CATS (1 offset, 0 no offset):
%           North  East   Up   Code
%             1      1     1     7   
%             1      1     0     6
%             1      0     1     5
%             1      0     0     4
%             0      1     1     3
%             0      1     0     2
%             0      0     1     1
%   (clearly, the triad 0-0-0, which corresponds to 0, is excluded)
% - If Mos is undefined or empty, no offsets are used.
%
% StrData is the (optional) input string for the kind of data to be passed
% to neu or mom files. Options:
% - if StrData is undefined, empty, is the string 'Raw' or is not an 
%   admitted string, the data are taken from the raw time series;
% - if StrData is 'Pre', the data are taken from the cleaned time series
%   (if no cleaned time series are unavailable, a warning message is shown,
%   no tasks are performed and tsOut = ts);
% - if StrData is 'PreCME', the data are taken from the preprocessed, 
%   CME-filtered time series (if these time series are unavailable, a 
%   warning message is shown, no tasks are performed and tsOut = ts).
%
% Mtsub is the (optional) matrix of time subdivision. Options: 
% - If Mtsub is a M-by-2 matrix Mtsub = [t1 t2], t1(k) and t2(k) are the 
%   minimum and maximum time limit respectively of the k-th partial time 
%   series in which the complete time series is subdivided. Please note 
%   that overlap between partial time series is allowed. These time limits 
%   are expressed either as MATLAB serial dates or dates in fractional year
%   form. 
%   Please note that a valid fractional data should be at least 1980.01366
%   and a valid serial date should be at least 723186 (6 January 1980);  
% - If Mtsub is a 2-elements vector, Mtsub(1) is the number M of output 
%   files. In this case, the data are subdivided into M partial time series
%   having equal length and whose overlap fraction is Mtsub(2). It is 
%   required that M is an integer and Mtsub(2) is in the range 0.1-0.5 (if 
%   M is not and integer, round(M) is used. If Mtsub(2) is lower than 0.1, 
%   0.1, i.e. 10% overlap, is used. If Mtsub(2) is higher than 0.5, 0.5,
%   i.e. 50% overlap, is used). A number of subdivisions such that each 
%   sub-time series is at least 3 years long is recommended. The initial 
%   date of the first sub-time series and the final date of the last 
%   sub-time series are the initial and final date of the whole time series 
%   respectively;  
% - If Mtsub is a scalar, it is the number M of output files. In this case,
%   the percentage of overlap is 30%. As in the previous case, round(Mtsub) 
%   is used;
% - If Mtsub is a 3-elements vectors, it is the number M of output files. 
%   In this case, the subdvision into M sub-time series with 30% overlap
%   is such that the initial date of the first sub-time series is Mtsub(2) 
%   and the final date of the last sub-time series is Mtsub(3). These dates
%   can be expressed either in MATLAB serial form or in year fractional 
%   form (please note that these dates could also be outside the limits of 
%   the whole time series);
% - If Mtsub is 4-elements vectors, Mtsub(1) is the number M of output 
%   files, Mtsub(2) is the overlap fraction (see above for the requirements
%   on this value) and Mtsub(3) and Mtsub(4) are the initial date of the 
%   first sub-time series and the final date of the last sub-time series 
%   respectively.
% - If Mtsub is undefined or empty, no a subdivision is carried out.
%
% The output logical msgOut is true if a subdivision is carried out, false 
% elsewhere.
% The output MtsubOut is the output M-by-2 matrix MtsubOut = [t1Out t2Out]
% with the really carried out subdivision (If Mtsub is a valid M-by-2
% matrix, it is MtsubOut = Mtsub).
%
% See also InspOffsetMom, tsData, tsDataFileIn.
%
%   [msgOut,MtsubOut,tsOut] = writeNeuMom(FILENA,OFP,tsIn,Mos,Mtsub)

% G. Teza, 2021, 2022.

if nargin < 6
    Mtsub = [];
end

if (nargin < 5) || isempty(strData) || ~ismember(strData,{'Raw','Pre','PreCME'})
    strData = 'Raw';
end

if nargin < 4
    Mos = [];
end

if nargin < 3 
    ts = [];
end

[ts,IT] = tsDataFileIn(ts); % tsData load

Is = searchFilledFields(ts);
if sum(Is) == 0 % all fields are empty - no neu/mom files can be generated
    IT = false;
end
msgOut = IT;

if ~IT
    MtsubOut = [];
    tsOut = [];
    return
end
tsOut = ts;

if (nargin < 2) || isempty(OFP) || (sum(logical(OFP)) == 0)
    OFP = false;
else
    if (numel(OFP) == 1) && (OFP == 2) 
        IEV = [0 0 1];  % scalar equal to 2: only vertical 
    elseif (numel(OFP) == 2) && (sum(OFP == [0 1]) == 2)
        IEV = [1 1 0];  % vector [0 1]: only horizontal
    else
        IEV = [1 1 1];  % all other cases: all components
    end
    OFP = true;
end

if nargin < 1
    filena = [];
end

% Information from tsData object:
statName = ts.statName;
t = ts.t;
dateyfrac = ts.dateyfrac;
se = ts.sE;          
sn = ts.sN;          
sv = ts.sV;
if strcmpi(strData,'Raw')
    ed = ts.Ed;
    nd = ts.Nd;
    vd = ts.Vd;
else
    if strcmpi(strData,'Pre') && ~isempty(ts.cleanedE)
        tn = ts.cleanedE.t;
        ed = ts.cleanedE.p;
        nd = ts.cleanedN.p;
        if ~isempty(ts.cleanedV)
            vd = ts.cleanedV.p;
        else
            vd = zeros(size(tn));
        end
    elseif strcmpi(strData,'PreCME') && ~isempty(ts.CMEfiltered)
        tn = ts.CMEfiltered.t;
        ed = ts.CMEfiltered.E;
        nd = ts.CMEfiltered.N;
        if ~isempty(ts.CMEfiltered.V)
            vd = ts.CMEfiltered.V;
        else
            vd = zeros(size(tn));
        end
    else
        tsOut = ts;
        msgOut = [];
        MtsubOut = [];
        warning('No %s time series available for station %s',strData,statName); 
        return
    end
    Iotn = ismember(t,tn);
    if sum(Iotn) < numel(t)
        t = t(Iotn);
        sn = sn(Iotn);
        se = se(Iotn);
        sv = sv(Iotn);
        dateyfrac = dateyfrac(Iotn);
    end
end
ti = t(1);
tf = t(end);
    
% data matrix:
Md = [dateyfrac nd ed vd sn se sv];    % note column order

if isempty(Mtsub)
    nfile = 1;
else
    if isscalar(Mtsub) || (isvector(Mtsub) && (Mtsub(1) < 1900))
        if isscalar(Mtsub)
            nfile = round(Mtsub);
            if nfile > 1
                frac = 0.3;
            end
        else
            nfile = round(Mtsub(1));
            if numel(Mtsub) == 2
                frac = Mtsub(2);
            elseif numel(Mtsub) == 3
                frac = 0.3;
                ti = Mtsub(2);
                tf = Mtsub(3);
            else
                frac = Mtsub(2);
                ti = Mtsub(3);
                tf = Mtsub(4);
            end
            if ti < 2100
                ti = frac2serial(ti);
                tf = frac2serial(tf);
            end
            if nfile > 1  
                if frac < 0.1
                    frac = 0.1;
                end
                if frac > 0.5
                    frac = 0.5;
                end
            end
        end
        if nfile > 1
            % It is L = (tf-t1)/[n-(n-1)*f],         ti |------------- ti+L
            % t2k = ti+L[k-(k-1)f], t1k = t_2k-L          ti+L-fL ------------- ti+2L-fL 
            L = (tf-ti)/(nfile-(nfile-1)*frac);     %            ti+2L-2fL -------------| tf=ti+3L-2fL    
            vt = (1:nfile)';    % such a vector must be defined to allow vectorial computations 
            t2 = (ti+L*(vt-(vt-1)*frac))';
            t1 = t2-L;
            Mtsub = [serial2frac(t1) serial2frac(t2)];  % serial2frac works on vectors!
        else
            Mtsub = [serial2frac(ti) serial2frac(tf)];  % for the possible limits
        end
    else 
        nfile = size(Mtsub,1);
        if max(Mtsub(:,1)) > 2100
            Mtsub = [serial2frac(Mtsub(:,1)) serial2frac(Mtsub(:,2))];
        end
    end
end

if isempty(filena)
    filenae = statName;
else
    filenae = filena;
end

% strC = {'N','E','U'};
strC = {'1','0','2'};   % 0: E (2nd col NEU), 1: N (first col NEU), 2: U (3rd col NEU)
if nfile == 1
    if ~OFP
        filenac = {[filenae '.neu']};
    else
        nfilee = 3;
        filenac = cell(nfilee,1);
        for kk = 1:3
            if IEV(kk)
                filenackk = [filenae '_' strC{kk} '.mom'];
                filenac(kk) = {filenackk};
            end
        end
    end
else
    if ~OFP
        nfilee = nfile;
        filenac = cell(nfilee,1);
        for k = 1:nfile
            filenac(k) = {[filenae '_' num2str(k) '.neu']};
        end
    else
        nfilee = nfile*3;
        filenac = cell(nfilee,1);
        for k = 1:nfile
            for kk = 1:3
                if IEV(kk)
                    filenac(k+kk-1) = {[filenae '_' num2str(k) '_' strC{kk} '.mom']};
                end
            end
        end
    end
end

for m = 1:nfile

    if ~isempty(Mtsub)
        t1m = Mtsub(m,1);
        t2m = Mtsub(m,2);
    else
        t1m = Md(1,1);
        t2m = Md(end,1);
    end
    
    Idm = (Md(:,1) >= t1m) & (Md(:,1) <= t2m);
    
    if sum(Idm) > 0
        
        % file open and header write:
        if ~OFP
            
            Mdm = Md(Idm,:);    % m-th submatrix
            filenam = filenac{m}; 
            fidm = fopen(filenam,'w');
            fprintf(fidm,'# Site : %s\n',statName);
        
        else
            
            if IEV(1)
                Mdm1 = Md(Idm,1:2);
                Mdm1(:,1) = frac2MJD(Mdm1(:,1));
                Mdm1(:,2) = 1000*Mdm1(:,2); % to have data in mm
                filenam1 = filenac{3*(m-1)+1};
                fidm1 = fopen(filenam1,'w');
                fprintf(fidm1,'# sampling period 1.0\n');
            end
            
            if IEV(2)
                Mdm2 = Md(Idm,[1 3]);
                Mdm2(:,1) = frac2MJD(Mdm2(:,1));
                Mdm2(:,2) = 1000*Mdm2(:,2); % to have data in mm
                filenam2 = filenac{3*(m-1)+2};
                fidm2 = fopen(filenam2,'w');
                fprintf(fidm2,'# sampling period 1.0\n');
            end

            if IEV(3)
                Mdm3 = Md(Idm,[1 4]);
                Mdm3(:,1) = frac2MJD(Mdm3(:,1));
                Mdm3(:,2) = 1000*Mdm3(:,2); % to have data in mm
                filenam3 = filenac{3*(m-1)+3};
                fidm3 = fopen(filenam3,'w');
                fprintf(fidm3,'# sampling period 1.0\n');
            end
            
        end

        if ~isempty(Mos)
            Iosm = (Mos(:,1) >= t1m) & (Mos(:,1) <= t2m);
            Mosm = Mos(Iosm,:);
            if ~isempty(Mosm)
                if ~OFP 
                    nosm = size(Mosm,1);
                    for k = 1:nosm
                        fprintf(fidm,'# offset %4.5f %d\n',Mosm(k,1),Mosm(k,2));
                    end
                else
                    if IEV(1)
                        Inosm1 = ismember(Mosm(:,2),4:7);
                        nosm1 = sum(Inosm1);
                        if nosm1 > 0
                            Mosm1 = Mosm(Inosm1,:);
                            Mosm1(:,1) = frac2MJD(Mosm1(:,1));
                            for k = 1:nosm1
                                fprintf(fidm1,'# offset %5.1f\n',Mosm1(k,1));
                            end
                        end
                    end
                    if IEV(2)
                        Inosm2 = ismember(Mosm(:,2),[2 3 6 7]);
                        nosm2 = sum(Inosm2);
                        if nosm2 > 0
                            Mosm2 = Mosm(Inosm2,:);
                            Mosm2(:,1) = frac2MJD(Mosm2(:,1));
                            for k = 1:nosm2
                                fprintf(fidm2,'# offset %5.1f\n',Mosm2(k,1));
                            end
                        end
                    end
                    if IEV(3)
                        Inosm3 = ismember(Mosm(:,2),[1 3 5 7]);
                        nosm3 = sum(Inosm3);
                        if nosm3 > 0
                            Mosm3 = Mosm(Inosm3,:);
                            Mosm3(:,1) = frac2MJD(Mosm3(:,1));
                            for k = 1:nosm3
                                fprintf(fidm3,'# offset %5.1f\n',Mosm3(k,1));
                            end
                        end
                    end
                end
            end
        end

        % data:
        if ~OFP
            fprintf(fidm,'%4.5f\t%f\t%f\t%f\t%f\t%f\t%f\n',Mdm'); % note "'"
            fclose(fidm);
        else
            if IEV(1)
                fprintf(fidm1,'%5.1f\t%f\n',Mdm1'); % note "'"
                fclose(fidm1);
            end
            if IEV(2)
                fprintf(fidm2,'%5.1f\t%f\n',Mdm2'); 
                fclose(fidm2);
            end
            if IEV(3)
                fprintf(fidm3,'%5.1f\t%f\n',Mdm3'); 
                fclose(fidm3);
            end
        end
        
    end
end

MtsubOut = Mtsub;

% tsOut upgrade (Note that MOS is NEU, not ENV!): 

if ~isempty(Mos)
    if ~OFP || (OFP && IEV(1))
        InosN = ismember(Mos(:,2),4:7);
        nosN = sum(InosN);
        if nosN > 0
            OffN = frac2serial(Mos(InosN,1));
            tsOut.offsetsN = OffN;
        end
    end
    if ~OFP || (OFP && IEV(2))
        InosE = ismember(Mos(:,2),[2 3 6 7]);
        nosE = sum(InosE);
        if nosE > 0
            OffE = frac2serial(Mos(InosN,1));
            tsOut.offsetsE = OffE;
        end
    end
    if ~OFP || (OFP && IEV(3))
        InosV = ismember(Mos(:,2),[1 3 5 7]);
        nosV = sum(InosV);
        if nosV > 0
            OffV = frac2serial(Mos(InosV,1));
            tsOut.offsetsV = OffV;
        end
    end
end